#########################################################################
### Kline, Bankert, Levitan, Kraft. 2017.                             ###
### "Personality and Prosocial Behavior: A Multilevel Meta-Analysis"  ###
### - model estimation                                                ###
#########################################################################

## load packages etc
rm(list=ls())
library(rstan)
library(dplyr)
library(haven)

## set options
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

## set condition to estimate all models
rerun <- T


### load and prepare datasets (wd is set in MASTER.R)

raw <- read_dta("mlma_data.dta") %>%
  mutate(svo2 = (svo2 - ifelse(svo2==1 & study %in% unique(study[svo1!=svo2]), .5, 0)) *
           ifelse(study %in% unique(study[svo1!=svo2]), 2, 1) # correct original svo2 var
         , svoD1 = as.numeric(svo1>=.5)             # convert DV to dummy (version 1)
         , svoD2 = as.numeric(svo1>.5))             # convert DV to dummy (version 2)

## data for all studies with personality data big5/HEXACO and MBTI (12 studies)
dat1 <- raw %>%
  select(study, extro, open, agree, conscien, payment, hilbig, benner, mbti, coop
         , svo1, svo2, svoD1, svoD2) %>%                          # select variables
  filter(study %in% 3:20) %>%                       # select studies
  mutate(study = as.numeric(as.factor(study))) %>%  # convert study no. to count
  na.omit()                                         # remove missings (in raw$study==19)

## data for big5/ HEXACO studies only (10 studies)
dat2 <- raw %>%
  select(study, extro, open, agree, conscien, neuro, payment, hilbig, benner, coop
         , svo1, svo2, svoD1, svoD2) %>%                          # select variables
  filter(study %in% c(3:8,12:20)) %>%               # select studies
  mutate(study = as.numeric(as.factor(study))) %>%  # convert study no. to count
  na.omit()                                         # remove missings (in raw$study==19)

## data for all studies, for incentivization-only analysis (15 studies)
dat3 <- raw %>%
  select(study, payment, hilbig, benner, coop
         , svo1, svo2, svoD1, svoD2) %>%                          # select variables
  mutate(study = as.numeric(as.factor(study)))      # convert study no. to count

## summary of datasets
summary(dat1)
summary(dat2)
summary(dat3)


### Models original submission

if(rerun == F) load("results.Rdata")

## Model 1: Big 5 + MBTI (12 studies)
if(rerun){
    X <- dat1 %>% select(extro, open, agree, conscien)
    Z <- dat1 %>% select(payment, hilbig, benner, mbti, coop) %>%
        aggregate(list(study=dat1$study), mean)
    Z <- cbind(1,Z[,-1])
    m1dl <- list(n = nrow(X), m = nrow(Z), nX = ncol(X), nZ = ncol(Z)
               , X = as.matrix(X), Z = as.matrix(Z)
               , svo = dat1$svo1, study = dat1$study)
    m1stan <- stan(file = "mlma.stan", data=m1dl, pars = c("beta","gamma"), iter=10000, thin = 10)
}
print(m1stan, pars=c("beta","gamma"))

## Model 2: Big 5 (10 studies)
if(rerun){
    X <- dat2 %>% select(extro, open, agree, conscien, neuro)
    Z <- dat2 %>% select(payment, hilbig, benner, coop) %>%
        aggregate(list(study=dat2$study), mean)
    Z <- cbind(1,Z[,-1])
    m2dl <- list(n = nrow(X), m = nrow(Z), nX = ncol(X), nZ = ncol(Z)
               , X = as.matrix(X), Z = as.matrix(Z)
               , svo = dat2$svo1, study = dat2$study)
    m2stan <- stan(file = "mlma.stan", data=m2dl, pars = c("beta","gamma"), iter=20000, thin = 20)
}
print(m2stan, pars=c("beta","gamma"))

## Model 3: Incentivization (15 studies)
if(rerun){
    X <- matrix(0, nrow=nrow(dat3))
    Z <- dat3 %>% select(payment, coop) %>%
    aggregate(list(study=dat3$study), mean)
    Z <- cbind(1,Z[,-1])
    m3dl <- list(n = nrow(X), m = nrow(Z), nX = ncol(X), nZ = ncol(Z)
               , X = as.matrix(X), Z = as.matrix(Z)
               , svo = dat3$svo1, study = dat3$study)
    m3stan <- stan(file = "mlma.stan", data=m3dl, pars = c("beta","gamma"), iter=10000, thin = 10)
}
print(m3stan, pars=c("gamma"))


### Robustness checks (focusing on model 1 for now)

## 1) varying study inclusion
if(rerun){
    m4stan <- list(NULL)
    for(i in unique(dat1$study)){
        dat4 <- dat1 %>% filter(study != i) %>%
            mutate(study = as.numeric(as.factor(study)))
        X <- dat4 %>% select(extro, open, agree, conscien)
        Z <- dat4 %>% select(payment, hilbig, benner, mbti, coop) %>%
            aggregate(list(study=dat4$study), mean)
        Z <- cbind(1,Z[,-1])
        m4dl <- list(n = nrow(X), m = nrow(Z), nX = ncol(X), nZ = ncol(Z)
                   , X = as.matrix(X), Z = as.matrix(Z)
                   , svo = dat4$svo1, study = dat4$study)
        m4stan[[i]] <- stan(file = "mlma.stan", data=m4dl, pars = c("beta","gamma")
                          , iter=10000, thin = 10)
    }
}
lapply(m4stan, print, pars=c("beta","gamma"))

## 2) random slope models
if(rerun){
    X <- dat1 %>% select(extro, open, agree, conscien)
    Z <- dat1 %>% select(payment, hilbig, benner, mbti, coop) %>%
        aggregate(list(study=dat1$study), mean)
    Z <- cbind(1,Z[,-1])
    m5dl <- list(n = nrow(X), m = nrow(Z), nX = ncol(X), nZ = ncol(Z)
               , X = as.matrix(X), Z = as.matrix(Z)
               , svo = dat1$svo1, study = dat1$study)
    m5stan <- stan(file = "mlma_slope.stan", data=m5dl
                 , iter=10000, thin = 10)
}
print(m5stan)


### Rescaling pro-social behavior / SVO variable (for original model 1-3)

## Model 1: Big 5 + MBTI (12 studies)
if(rerun){
    X <- dat1 %>% select(extro, open, agree, conscien)
    Z <- dat1 %>% select(payment, hilbig, benner, mbti, coop) %>%
        aggregate(list(study=dat1$study), mean)
    Z <- cbind(1,Z[,-1])
    m1svo2dl <- list(n = nrow(X), m = nrow(Z), nX = ncol(X), nZ = ncol(Z)
                   , X = as.matrix(X), Z = as.matrix(Z)
                   , svo = dat1$svo2, study = dat1$study)
    m1svo2stan <- stan(file = "mlma.stan", data=m1svo2dl, pars = c("beta","gamma")
                     , iter=10000, thin = 10)
}
print(m1svo2stan, pars=c("beta","gamma"))

## Model 2: Big 5 (10 studies)
if(rerun){
    X <- dat2 %>% select(extro, open, agree, conscien, neuro)
    Z <- dat2 %>% select(payment, hilbig, benner, coop) %>%
        aggregate(list(study=dat2$study), mean)
    Z <- cbind(1,Z[,-1])
    m2svo2dl <- list(n = nrow(X), m = nrow(Z), nX = ncol(X), nZ = ncol(Z)
                   , X = as.matrix(X), Z = as.matrix(Z)
                   , svo = dat2$svo2, study = dat2$study)
    m2svo2stan <- stan(file = "mlma.stan", data=m2svo2dl, pars = c("beta","gamma")
                     , iter=20000, thin = 20)
}
print(m2svo2stan, pars=c("beta","gamma"))

## Model 3: Incentivization (15 studies)
if(rerun){
    X <- matrix(0, nrow=nrow(dat3))
    Z <- dat3 %>% select(payment, coop) %>%
    aggregate(list(study=dat3$study), mean)
    Z <- cbind(1,Z[,-1])
    m3svo2dl <- list(n = nrow(X), m = nrow(Z), nX = ncol(X), nZ = ncol(Z)
                   , X = as.matrix(X), Z = as.matrix(Z)
                   , svo = dat3$svo2, study = dat3$study)
    m3svo2stan <- stan(file = "mlma.stan", data=m3svo2dl, pars = c("beta","gamma")
                     , iter=10000, thin = 10)
}
print(m3svo2stan, pars=c("beta","gamma"))


### Rescaling pro-social behavior / SVO variable to dummy (for original model 1-3)

## Model 1: Big 5 + MBTI (12 studies)
if(rerun){
    X <- dat1 %>% select(extro, open, agree, conscien)
    Z <- dat1 %>% select(payment, hilbig, benner, mbti, coop) %>%
        aggregate(list(study=dat1$study), mean)
    Z <- cbind(1,Z[,-1])
    m1svoD1dl <- list(n = nrow(X), m = nrow(Z), nX = ncol(X), nZ = ncol(Z)
                   , X = as.matrix(X), Z = as.matrix(Z)
                   , svo = dat1$svoD1, study = dat1$study)
    m1svoD1stan <- stan(file = "mlma_logit.stan", data=m1svoD1dl, pars = c("beta","gamma")
                     , iter=10000, thin = 10)
    m1svoD2dl <- list(n = nrow(X), m = nrow(Z), nX = ncol(X), nZ = ncol(Z)
                   , X = as.matrix(X), Z = as.matrix(Z)
                   , svo = dat1$svoD2, study = dat1$study)
    m1svoD2stan <- stan(file = "mlma_logit.stan", data=m1svoD2dl, pars = c("beta","gamma")
                     , iter=10000, thin = 10)
}
print(m1svoD1stan, pars=c("beta","gamma"))
print(m1svoD2stan, pars=c("beta","gamma"))

## Model 2: Big 5 (10 studies)
if(rerun){
    X <- dat2 %>% select(extro, open, agree, conscien, neuro)
    Z <- dat2 %>% select(payment, hilbig, benner, coop) %>%
        aggregate(list(study=dat2$study), mean)
    Z <- cbind(1,Z[,-1])
    m2svoD1dl <- list(n = nrow(X), m = nrow(Z), nX = ncol(X), nZ = ncol(Z)
                   , X = as.matrix(X), Z = as.matrix(Z)
                   , svo = dat2$svoD1, study = dat2$study)
    m2svoD1stan <- stan(file = "mlma_logit.stan", data=m2svoD1dl, pars = c("beta","gamma")
                     , iter=20000, thin = 20)
    m2svoD2dl <- list(n = nrow(X), m = nrow(Z), nX = ncol(X), nZ = ncol(Z)
                   , X = as.matrix(X), Z = as.matrix(Z)
                   , svo = dat2$svoD2, study = dat2$study)
    m2svoD2stan <- stan(file = "mlma_logit.stan", data=m2svoD2dl, pars = c("beta","gamma")
                     , iter=20000, thin = 20)
}
print(m2svoD1stan, pars=c("beta","gamma"))
print(m2svoD2stan, pars=c("beta","gamma"))

## Model 3: Incentivization (15 studies)
if(rerun){
    X <- matrix(0, nrow=nrow(dat3))
    Z <- dat3 %>% select(payment, coop) %>%
    aggregate(list(study=dat3$study), mean)
    Z <- cbind(1,Z[,-1])
    m3svoD1dl <- list(n = nrow(X), m = nrow(Z), nX = ncol(X), nZ = ncol(Z)
                   , X = as.matrix(X), Z = as.matrix(Z)
                   , svo = dat3$svoD1, study = dat3$study)
    m3svoD1stan <- stan(file = "mlma_logit.stan", data=m3svoD1dl, pars = c("beta","gamma")
                     , iter=10000, thin = 10)
    m3svoD2dl <- list(n = nrow(X), m = nrow(Z), nX = ncol(X), nZ = ncol(Z)
                   , X = as.matrix(X), Z = as.matrix(Z)
                   , svo = dat3$svoD2, study = dat3$study)
    m3svoD2stan <- stan(file = "mlma_logit.stan", data=m3svoD2dl, pars = c("beta","gamma")
                      , iter=10000, thin = 10)
}
print(m3svoD1stan, pars=c("beta","gamma"))
print(m3svoD2stan, pars=c("beta","gamma"))


### Robustness check: drop cooperation dummy in Model 1-3

## Model 1: Big 5 + MBTI (12 studies)
if(rerun){
    X <- dat1 %>% select(extro, open, agree, conscien)
    Z <- dat1 %>% select(payment, hilbig, benner, mbti) %>%
        aggregate(list(study=dat1$study), mean)
    Z <- cbind(1,Z[,-1])
    m1coop0dl <- list(n = nrow(X), m = nrow(Z), nX = ncol(X), nZ = ncol(Z)
               , X = as.matrix(X), Z = as.matrix(Z)
               , svo = dat1$svo1, study = dat1$study)
    m1coop0stan <- stan(file = "mlma.stan", data=m1coop0dl, pars = c("beta","gamma")
                      , iter=10000, thin = 10)
}
print(m1coop0stan, pars=c("beta","gamma"))

## Model 2: Big 5 (10 studies)
if(rerun){
    X <- dat2 %>% select(extro, open, agree, conscien, neuro)
    Z <- dat2 %>% select(payment, hilbig, benner) %>%
        aggregate(list(study=dat2$study), mean)
    Z <- cbind(1,Z[,-1])
    m2coop0dl <- list(n = nrow(X), m = nrow(Z), nX = ncol(X), nZ = ncol(Z)
               , X = as.matrix(X), Z = as.matrix(Z)
               , svo = dat2$svo1, study = dat2$study)
    m2coop0stan <- stan(file = "mlma.stan", data=m2coop0dl, pars = c("beta","gamma")
                 , iter=20000, thin = 20)
}
print(m2coop0stan, pars=c("beta","gamma"))

## Model 3: Incentivization (15 studies)
if(rerun){
    X <- matrix(0, nrow=nrow(dat3))
    Z <- dat3 %>% select(payment) %>%
    aggregate(list(study=dat3$study), mean)
    Z <- cbind(1,Z[,-1])
    m3coop0dl <- list(n = nrow(X), m = nrow(Z), nX = ncol(X), nZ = ncol(Z)
               , X = as.matrix(X), Z = as.matrix(Z)
               , svo = dat3$svo1, study = dat3$study)
    m3coop0stan <- stan(file = "mlma.stan", data=m3coop0dl, pars = c("beta","gamma")
                      , iter=10000, thin = 10)
}
print(m3coop0stan, pars=c("gamma"))




### save model results

if(rerun) {
    save(list = ls()[c(grep("stan",ls()), grep("dat",ls()))], file = "results.Rdata")
}


